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ABSTRACT 


This report presents a new technique for the design of 
multidimensional FIR digital filters. Given the desired amplitude 
response, the truncation error of a low order approximation to a high 
order filter is iteratively decreased by varying the phase response. 
Convergence of the iterative algorithm is discussed. In general, the 
phase responses of the filters resulting from the algorithm are nonlinear. 
Pictures suffer phase distortion when filters with nonlinear phase char- 
acteristics are used in filtering. A measure of phase distortion in 
pictures is introduced so that the effects of nonlinear phase in filters 
may be quantified. Examples illustrating the use of the design tech- 
nique and the measure of phase distortion are included. The results 
of this report are formulated in two dimensions; however they are 


valid for filters and signals of any dimension. 
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I. INTRODUCTION 


A. Two-Dimensional Digital Filtering 


In the last several years digital techniques have found 
widespread use in image processing due to the availability of digital 
computers, the large number of operations which computers can per- 
form, the increasing speed at which these operations can be performed, 
and the advent of efficient algorithms for carrying out certain process- 
ing tasks. The operation most widely used in digital image processing 
is the linear operation known as convolution or filtering. In [1], Rind- 
fleisch, Dunne, Frieden, Stromberg, and Ruiz describe the use of con- 
volution and other digital processing techniques for the enhancement 
of the Mariner 6 and 7 pictures of Mars. Many degradations of pic- 
tures, such as motion and defocus blurs, may be described as the 
convolution of an ideal picture with a degrading signal. Stockham, 
Cannon, and Ingebretsen [2] have shown how to identify the degrading 
signal, when it is not known, and how to remove it through a second 
convolution, which is known as a deconvolution. Lewis and Sakrison 
[3] have used Wiener filters to reduce noise and optical distortion in 
electron micrographs. Duda and Hart [4] have suggested filtering as 
an aid to scene analysis. They have shown the effects, on pictures, 
of low pass, high pass, and high emphasis two-dimensional digital 
filtering. In the oil industry two dimensional digital filtering is being 
used to remove multiple reverberations from marine seismic data. 

Let a be a digital input sequence or picture and let 


fh? be a two-dimensional impulse response sequence. We assume 


that fh! has M rows and N columns and that han! is the same 


size or larger. When the input sequence ae. is convolved with 

or filtered by the sequence i) , the output picture (Bn! is 

described by * 
a M-1 N-1 


= z= . 
Pan x20 220 "xe! ree) 


kt m-k,n-t 
The two-dimensional digital filter represented by ih) may also be 
described by the Z-transform transfer function, 


M-1 N-1 


Ha (Z, .Z,) = (152) 


Ha (2, .2,) is calied a finite impulse response or FIR digital filter if 
M and N are finite nonnegative integers. If M or N are infinite 
positive integers, Ha (2) +2.) is called an infinite impulse response 
or IIR digital filter. The frequency response of the filter Ha (2, 2) 
is Ha (exp(-j#,), exp(-j,) ), which may be evaluated through use 
of the discrete Fourier transform or DFT. For every spatial sequence 
fh! of length (M,N), we can find a corresponding frequency do- 
main sequence {H,,}, which is the DFT of ae The DFT rela- 
tionship, denoted as hon Fy ys is described as 


M-1 N-1 k 
A= 3 8. si a (1.3) 
m=0 n=0 
and 
M- - a “ 
n= a/MN) MS? Nsw uty, aay 
k=0 g-0 “4 


where U = exp(-j 2™/M) and V=exp(-j 27/N). Thus a filter des- 
cribed by (1.2) is specified by the MN points of the spatial sequence, 


the, }, or by the MN points of the frequency domain sequence, {H,}. 


The frequency domain sequence is expressed by an amplitude sequence, 


(ln, 13, and a phase sequence, {6 3, where 
Se : 1S 
Hy = [Hy yl expe, (1.5) 
In general, for a given amplitude sequence, distinct choices of the 


phase sequence {e,,}, lead to distinct impulse response sequences, 


{hat The relationship between Ha (Z,.Z,) and {HL} is 


ee ~ He (exp(-jw,),exp(-jw,) ) (1.6) 


where 


(w +@,) = (@T/M)k, (27/N)2) . (1.7) 


B. FIR Digital Filter Design 


In designing a two-dimensional digital filter we are 
usually given a desired frequency response which may be expressed, 


using amplitude and phase, as 


- gd 
H(w,,W,) = |H(w,,w,)lexp Ge (w),w,)) . (1.8) 
One may design an FIR digital filter using (1.4), (1.7), and 


H = Hw, ,)) : (1.9) 


kt 


The resulting filter Ha (2, 2) has the desired frequency response at 


the MN points of (1.7). However, if M and N are large enough for 


Ha (exp(-j,), exp (-j#,)) to closely approximate H(w,,,), then 

Ha (2, +2) usually has too many coefficients to be useful. This 

problem is attacked by designing a filter Ha (2) 4%), with few coef- “ 
ficients, which approximates Ha (2, +2) or H(w),,). Let Ty be 

the set of nonnegative integer pairs (m,n), where 0 =m* M-1l and 


0 =n=N-1, and let T bea subset of T.. Given the spatial se- 


0 
quence {he (m,n) eT). we wish to form a filter 
= -l*s m on 
mee oye ORE eee (1.10) 
aye mn ~1 Z 
m=0 n=0 


such that ee =0 if (m,n) ¢T. The other coefficients are chosen so 
that the approximation is accurate. In general T may be any set of 
nonnegative integer pairs satisfying the condition above, however in 


this report the choice of T is restricted to 


T= {(m,n) ¢2 |0<msM,,0<n=N,}, (eld) 


1 


where M, <M-l, Ny <N-1l, and Z denotes the set of all integer 
pairs. In this case Ha (2, +2.) is a low order approximation of 


Ha (2) 1Z,). 


Design Techniques 
Recently minimax methods for designing two-dimensional 
FIR filters have been explored by Hersey and Mersereau [5]. Their 


technique minimizes the error function 


En i a ae oy - | Ha (exp(-j4,) ,exp(-ja,))|)| . 


where Ww, .,) is a weighting function and Ha (2, +25) is restricted 


to linear phase. Unfortunately minimax design techniques require a 
prohibitive amount of computation time. 

Another two-dimensional FIR filter design method is the 
frequency sampling method developed by Hu and Rabiner [6]. When 
a filter is designed by frequency sampling the DFT frequency domain 
coetficients, corresponding to pass-bands and stop-bands, are fixed. 
The error of approximation is minimized with respect to the transition- 
band DFT coefficients. The impulse response coefficients are then 
calculated by taking the inverse DFT. An advantage of this filter de- 
sign method is that it can be used to approximate any frequency re- 
sponse. This method consumes a lot of computation time, but not as 
much time as a minimax method. 

Two-dimensional FIR digital filters may be designed using 
the transformation approach of McClellan [7]. In this approach the 


equation 


cos (w) = Acos(#,) + Bcos(w,) - C cos(,)cos(w,) +D 


is used to transform a one-dimensional filter H (w) into a two- 
dimensional filter Hw, w,) . Different types of filters may be de- 
signed by changing the constants A, B, C, and D. This is the fast~ 
est of the design techniques. However, the method applies only to 
linear phase filters and there are severe amplitude constraints as well. 
Truncation is the simplest design method for two- 
dimensional FIR digital filters. Given Ha (2, .2,) and the set T, 
Ha (2, +2,) is formed using 


h  =h ,(m,n)eT , 
mn mn 


=0 ,(m,n)%T . 


Let {H,,} be the order (M,N) DFT of he. Truncation minimizes 
the summation of squares error function, 
M-1 N-1 


= ee 
Ee ey ey ae Bal 


; (1.13) 
with respect to oe This filter design method has the advantage 
that any desired frequency response may be approximated. There are 
no amplitude or linear phase constraints. Also, truncation is an ex- 
tremely fast design technique. Unfortunately the frequency domain 
error, due to truncation, is usually severe near transition regions. 
To lessen the bad effects of truncation, windowing is 

often employed. When (1.12) is replaced by 


h = 


h_ow__,(m,n)eT , 
mn mn mn 


(1.14) ‘ 
=0 , (m,n)sT , 
where {wind is a window sequence, the severe ripples in the fre- t 
quency response of the truncated filter are smoothed. Huang [8] has 
demonstrated that two-dimensional windows may be designed from one- 
dimensional windows. A problem with windowing is that it extends 
transition regions. This is undesirable when sharp cutoffs are re- 
quired in the amplitude response. In addition windowing is not an 
optimal design method. 
Of all the two-dimensional FIR digital filter design 
methods, only truncation and windowing lead to filters, which ap- 
proximate any desired frequency response, in a reasonable amount 
of computation time. The field of image processing may be signifi- 


cantly advanced if the truncation design method is improved and the 


need for windowing is reduced or eliminated. 


sida cli ri lan 


C. Phase Distortion 


It is being increasingly recognized that when complex 
signals such as Fourier transforms are represented in polar coor- 
dinates, phase is more important than amplitude. Larmore, El-Sum, 
and Metherell [9] and also Kermisch [10] have shown the importance 
of phase information in holography. Huang, Burnett, and Deczky [11] 
have demonstrated the significance of phase in digital image process- 
ing filters. In general the phase of a digital filter will have error. 
Linear phase, which does not degrade the phase of input signals, 
cannot be obtained in causal, realizable IIR digital filters. Even 
though FIR digital filters may have linear phase, this phase often 
has discontinuities of ™. Nonlinear phase characteristics, which 
may be specified for Wiener filters and matched filters, are gener- 
ally not realizable with zero error. The filter design method of 
Steiglitz [12] minimizes a magnitude error function and leads to 
nonlinear phase characteristics. 

let H(w, ,W,) be the frequency response of a digital 
filter. Let H(w).,) be the desired frequency response, as in (1.8). 
H(w,,¥,) has amplitude and phase error and may be the frequency 
response of an IIR or FIR digital filter. The frequency response may 


be expressed, using amplitude and phase, as 
a - a 8 
H(w,,w,) = |H(w,,w,)lexp(je(w,,,)) . (1.15) 


The Fourier transform of a digital input picture Uta! with M rows 


and N columns is 


M-1 N-1 


uf f__exp(-j(mw 
m=0 n=0 mn xP( i( 1 


F(w, w) = +nw,)) . (1.16) 


If fe} is filtered using the desired filter the output picture {pnt 
has the Fourier transform 
P(w, ,W,) = F(W,, W))H(w, 1 W,) : 
(1.17) 
By d 
= F(w, , @,) | H(w, , w,) lexp (je (w.@,)) 


The effects of the phase error of H(w, ,W,) are of interest to us, so 
we therefore construct a frequency response 

H(w,,W,) = [Hw ,w,)lexpG2(w,,,)) , (1.18) 
which has the phase error but not the amplitude error of H(w, ,,) : 


When H(w,,@,) is used to filter the input picture it te the out- 


put picture is Se which has the Fourier transform 


. #3 = ig 
P(w, ,W,) F(W), @,)H(, ,W,) 
= P(w) ,W,)expm(w, ,w,)) (1.19) ° 
where 
ye _ ae | 
m(W,,W,) = 8(w, ,w,) 8 (w,-W,) ; (1.20) 


The error between {pont and te! is called phase distortion and 

is caused by the phase error function m(W,,W,). The phase response 

8(w, ,W,) is acceptable in a given application if the phase distortion 

produced by m(W, ,W,) is measured and found to be small. One may 

attempt to measure the phase distortion of filtering by measuring 

m(w, 1 Wy) - However, such a measure does not change as the input 

picture changes. This type of measure has the same value for sine , 


wave and wideband input pictures. 
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D. Objectives 


The objectives of this report are 

(1) To develop a two-dimensional FIR digital filter 
design method which leads to amplitude responses, 
better than those obtainable through truncation and 
windowing, 

(2) To develop a measure of phase distortion, caused 

by filtering, for the evaluation of the phase responses 
of filters. 

A design method for two-dimensional FIR digital filters, 
which we call phase correction, is developed and analyzed in Chap- 
ter II of this document. The basic design algorithm is derived and 
its convergence is discussed in section IIA. In section IIB the algo- 
rithm is improved so that convergence is more rapid. A comparison of 
phase correction to gradient techniques is made in section IIC. The 
problems of initializing the algorithm and measuring the phase error of 
the filters are treated in section IID. 

In Chapter III a measure of phase distortion in pictures is 
formulated and its uses are described. The measure of phase distor- 
tion is developed from the Euclidean distance between two pictures in 
section IIIA. Methods of calculating the measure and evaluating phase 
responses of filters are given in section IIIB. An efficient way to es- 
timate the measure is formulated in section IIIC. 

We design several filters in section IVA to illustrate the 
use of the phase correction algorithm. Filters designed by phase cor- 
rection are compared to filters designed by truncation and windowing. 
In section IVB two filters are used to fiiter a test pattern input picture 
and the resulting phase distortion in the output pictures is measured 


and estimated. 
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II. THE PHASE CORRECTION ALGORITHM 


A. Derivation of the Algorithm 


Given the desired amplitude sequence fla, and the 
set T, the truncation error between Ha (2, +2) and Ha (2, .2,) may be 


written as 
M-1 N-l 
iia. = = ih 
m=0 n=0 
(m,n) £T 


mn| 


2,.2,. M-1 N-] ,M-1 N-1 = - 
=—Q/QAN ) f£ “xr [re = [H, ,lexpUe, ,)U my a os 
m=0 n=0 k=0 £0 


(m,n) £T 
a function of the phase sequence fe 3 . The expression (2.1) is 


(2.1) 


arrived at by the use of the earlier equations (1.2), (1.4), (1.5), ‘ 
and (1.10). (18, 3) is a measure of the goodness of the approxi- 
mation of Ha (2, +2) by Ha (2, .2,). In this section we derive an 
iterative method for picking phase sequences, {@, 13 , so that the 
error function 1({8, , }) is decreased with respect to {6 3 : 
For the ith iteration, where i is a nonnegative integer, the 
high order filter and its truncated version are described respectively 


as 


i a a 
Ha (2, ,Z,) = “shoe 2 


(2.2) 


11 


where 


ai i ; 
h = ian! (m,n) el, 


=0  , (m,n) €T. (2.4) 


It may be emphasized that i is a superscript and not an exponent. 
The frequency domain sequences, (H,,} and (a, pie are obtained 


from spatial sequences as 


i i 2.5) 
ON aaa te r ( 
and 
ai ai 
mn kg 


The amplitude and phase relationships are 


i i 
Big IHL | expe, ) (2.7) 
and 
at = (a! | expé ) (2.8) 
ee Ae, kp? - 


As seen in (2.7) the amplitude sequences (lay, 1 are restricted by 
the equation ln, = | for every 1, which implies that the am- 
plitude sequence is held to be the same for every iteration. Further, 


a simple application of Parseval's equation leads to 


M-1 N-l M-1 N-1 


1 2 
Ee th) = MN) METNST i y?ez (2.9) 
m=0 n=0 k=0 go “4 


12 


and 
nt ys ee (2.10) 
mn 
for every i, m, and n. 
Four equivalent expressions for the error between HA (2, Z,) 


and He (2) ,2,), where i and p are nonnegative integers, are: 


TT 7 
(i) r= (1/4n”) J J | Halexp(-tw,), exp(-jw) ) 


2 
- Ha (exp(-jo,) ,exp(-Jwy) ) |” dw, dw, (2.11) 
GS ee ee (2.12) 
mn mn 
m=0 n=0 
ait) TP = ymin) ME? N5} Le - Le (2.13) 
k=0 g¢0 
and 
i M-1 N-1 a ae 
(iv) PP = mn) ME} NS + TERI 
k=0 g¢0 
r7P » . 9 
-2 1H | [HY ,| cos(a, | a) (2.14) 


Expressions (2.12) and (2.14) follow directly from the evaluation of 
(2.11) and (2.13), respectively. Expression (2.13) follows when 


i -* He i ie 
Sn MN SS (2.15) 
k=0 0 


and 


13 
“ M-1 N-1 + ee 
he = (/MN) s x AP ou km gn ae 
i ia L 

k=0 g¢0 


are substituted into (2.12). Using (2.4) in (2.12) and letting p =i, 
we find that 


fi M-1 N-1 ,1i 2 
r=: = ho .) (2.17) 
mn 
m=0 n=0 
(m,n) £T 


and using (2.1), (2.5), (2.7), and (2.17), that 
Ufo J) aa (2.18) 


From (2.18) it is obvious that an algorithm which iteratively produces 


@ decreasing error sequence, 


also decreases the function 1(£8 , 3) iteratively. 

Given the desired amplitude sequence (lH, ,| } and the 
set T, the first step of the phase correction algorithm is to choose an 
initial phase sequence Lhe In order for the sequences (hr) and 
fh! to be real, all phase sequences, including the initial phase 
sequence, must satisfy the symmetry condition, 


i 


2 ke &f 2. 
y= FMM, 9(M-k]  /IN- 2], ’ (2.19) 


where ny.) is an integer sequence, and [M-k],,. [N-t] are in- 
iv 


tegers modulo M and N respectively. 
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The numbers I, and Ae are specified, in the algorithm's 


first step, to provide an aon from the algorithm. The integer eee is 
the maximum value of i allowed. In other words, only (i ax +1) iter- 
ations are allowed. - is a nonnegative real number such that if 
ph - 7p < I, for some i, the algorithm is halted. The steps of 


the algorithm are given below. 


Step 0: Given {Hy 1}. (ey 5} T, Ine and 4 ax? compute (Hy) 
from (2.7). 

Step 1: Set i equal to 0. 

Step 2: Calculate ihe using the relationship (2.5). 

Step 3: Get (hr from (2.4). 

Step 4: Calculate Ne using (2.17). It may be noted, from (2.12) 
and (2.17), that 


os epi} (2.20) 


i-l,i-1 Jil oe 
-I <I,ort=i) ; 


Exit if I ane 


Step 5: Compute (8 5] from (2.6) and (2.8). 
Step 6: Increase the index i by l. 
Step 7: Using the equation 
i i-l 
calculate (H,,) from (2.7). 


Step 8: Goto Step2. 


i-l 


i-l,i-1 
can be expressed by (2.14), pe can 


Whereas I 
be expressed, using (2.14) and (2.21), by 


15 


zs 1 N=1 ed-1 4,2 2.22 
fel _ ayn) ME (Hy yt - Fy )- aon 
k=0 ei 
This means that 
perl < poled), (2.23) 


The implementation of the phase correction algorithm is 
shown in Fig. 1. In the figure we have implemented Step 7 with 


Phas ai-1 , Ai-1 
H, = 1H, @,, 7 By, I) (2.24) 


which follows from (2.8), since 


gil) 8 -F ; 
expe, ,) = re Lan : (2.25) 


By alternately decreasing the truncation error in the 
spatial and frequency domains, the algorithm generates new spatial 
and frequency domain sequences for each iteration. Inequalities 
(2.20) and (2.23) imply that the sequence P is nonincreasing, so 
we have shown that the phase correction algorithm decreases 1(f@) )) 
or leaves it the same. 


If i =o and I_ =0 the infinite sequence 
max E 


00 1 
We ef et (2.26) 


is generated. 


+ + +1,i+ 
THEOREM 1: If rf = ? - or ? eine 2 naete for some finite 


nonegative integer i, the sequence W converges ina finite number 


of terms. 


1 


ii + 
PROOF: Case fia thet 


By using (2.13) and (2.22) one obtains 

2 
M=1 N=1 \at By amet Se (UH, l- eee (2.27) 
k=0 pit k=0 0 


; 2 
This implies that (or) = (e3 and that ly — {et }. Since the 
same operations are performed in each iteration, (oP) = (er, } for 


every nonnegative integer p. 


Case 2 pe = geet. 


By using (2.22) one obtains 


1 
pri.d (iymn) MS N- =) | lH, | expe) 
k=0 g=0 

il i+], 2 , 
- [H, | exp Ge I 
, Mad Netz del ot 2 
ae z Le se 

m9 n=0 


M=1 N-1 ,i+1 i \2. M=1 N-1 ,i+1 2 
Se ee = + 


Red = S hh y (2.28) 
m=0 n=0 as = m=0 n=0 sis 
(m,n) ¢T (m,n) #T 


Now Phe PAE ig .17) and @.26) imply that the right hand 


side of (2.28) equals 
M-1 N-1 i+1 2 


b> Se te) 
m=0 n=o [™ 
(m,n) £T 
and that 
M-l N-1,i1 «i (2 
r z Vin “ha =i, (2.29) 
m=0 n=0 


(m,n) ¢T 


| 
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Equation (2.29) implies that (hd = (A, af and (a7) = (é gis 
+ 
which leads to ter) = tet y. The same OE are Pits out 


kt 

+14 +1 
for each iteration un therefore ic, atp j= (ery } for every nonnega- 
tive integer p. 


Q.E.D. 


Corollary: The sequence W is always convergent. 
Proof: This is proved for a special case in theorem 1. From [13] all 
bounded, non-increasing sequences of real numbers converge. The 
members of W are bounded by 0 and and W is non-increasing 
as seen in (2.20) and (2.23). The corollary follows. 

Q.5.D. 


The convergence of the error sequence W is proved above, 


but this does not imply the convergence of the sequences, 


Sok a). eens | ee 
and 
Ye{ (h.), [Oe ey oak Eee 


The convergence properties of X and Y will be discussed ina paper. 


B. Increasing the Rate of Convergence 


Normally the phase correction algorithm reduces the error 


function 1({6, 2) significantly in the first few iterations. However, 


there is often room for improvement. In this section a modified 
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algorithm, which increases the rate of convergence of the error 
sequence, is presented. In particular, an algorithm that approxi- 


mates the high order filter 


fia, mies BSG 2s (2.30) 
2 mn 1 2 
m=0 n=0 
by the truncated filter 
-l1 N-1l-i mn 
Hh (2 a= ™%% sh. 2. 2, » (2.31) 
{2 mn 2 
m=0 n=0 
where the conditions 
~-i i 
= A (m,n) eT 
= 0, (m,n) £T @.32) 
and 
fi esin expt.) 
D inn | kyl P Fey (2.33) 


are satisfied, is described in the ee The phase sequences of 
the high order filters are denoted by (4) as in (2.33). As with the 


original algorithm every high order filter Ai (2, 2, ) has {lH,_,|} as 


kt, 
its amplitude sequence. An expression for a error between Ha (2) se) 


2 
and HB (2, -Z,) is 
~ip_ M-1 N-l i =p 2 
ic _ 2.34 
I bs : (Bon Ba * ( ) 


m=0 n=0 


The error sequence of the modified phase correction algorithm is given 
by 


SEE vasksl Et sawed (2.35) 


= 


In a manner analogous to the original algorithm, the first spatial 


sequence is generated from {lH 13 and {oP}. using the equation 


(i |) = IDFT (1H, JexpGg? )}. (2.36) 


Further, the first two iterations of the modified algorithm are identical 


to the first two iterations of the algorithm of section IIA, meaning that 


for i=0,1. For i=2 the equations above also hold, but for i22 we 


i 
form a second truncated sequence ee in each iteration by using 


=] 
mn 


. 3 Se ok 
n mn 


” ) (2.37) 


where p! is a real number chosen to accelerate the convergence of the 
algorithm. It may be assumed that in the first steps of the present al- 
gorithm, bp” and b are zero. The number pt is chosen so that the 
sum of squared errors, 


i 6 


f in? 2 
mn (h 


= M21 N-1 
nfo ae 


m=0 n=0 


= i-] 1—-i-l - 
K = 
i (h mn* b Mi ans h 


between {hi } and fit? + ptqin? . B43 is minimized. 
mn mn mn mn 
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achieved by differentiating (2.38) with respect to bi and setting the 


result equal to zero. This gives 


M-l Nel <i. ci-l,-i- - 1-2 
| a z Gon Flan Bae 
i i m=0 n=0 
tig ce oe oe 2.2 ; gh 
j - -1 ,<i-l -i- 
| py E &- Aan) 
i m=0 n=0 
Given the sequence ty’ the operations 
oe i (2 .40) 
{H,. a a Rs Sere cs | 
] 
ee a (2.41) ] 
~ —i-] = i-1 
AT |H, | (Hy / re (2 ,42) | 
and 
fl 3 = IDFT ea } (2.43) ‘ | 
{ mn kg 4 
are carried out for the implementation of the algorithm. This leaves the 
algorithm at the truncation step. Because of the way the sequences 
ty) are generated one cannot guarantee that 
pil er} a (2.44) 


Therefore, if (2.44) doesn't hold, let 


it «fi, (2.45) 


B= 0, (2.46) 
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and continue the iteration. In other words, if an iteration of the modi- 
fied algorithm fails to produce a smaller error number than the previous 
iteration, we repeat the iteration with . = 0, which is equivalent to 
performing an iteration in the original algorithm. The final error se- 
quence Ww has the property that 

eee rit (2.47) 
for all i greater than 0. If the iteration fails in the above fashion the 
value of the error function for two successive iterations is the same, 


meaning 


oe ee (2.48) 


However, this does not imply that the modified algorithm has con- 
verged. Whenever (2.45) is implemented, the equation 

={-] <-i-2 

h mn : mn 
holds in the following iteration, giving a zero denominator in (2.39). 
This is easily remedied by replacing ae by i in (2.39), in an 
iteration following the implementation of (2.45). 

The block diagram for the modified phase correction algo- 

rithm is shown in Figure 2. As with the original algorithm, the num- 
bers I, and lag provide exits. The steps of the algorithm are given 


below in detail. As seen in (2.39), be can't be calculated if i <2. 


Therefore, (2.4€‘ is used for i=0 and l. 


ee 
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. 0 
Step 0: Given lH, 41} fp} T, In, andi, set i equal to 


‘, 


zero and compute (a? i using 


| 
| 


Ne 


? 0 
Hy y= [H, ,lexpde, )- (2.49) 


Step 1: Calculate (fi using (2.43). 
Step 2: If (2.44) holds, continue to Step 3; otherwise skip to Step 6. 


Step 3: Perform the truncation in (2.32). 


pi-li-1 _ it 


Step 4: Exit if I <I, ori2i- ‘ 


ax 
Step 5: Ifi <2 set b! = 0; otherwise calculate b! from (2.39); go 

to Step 7. | 
Step 6: Perform (2.45) and (2.46). 


Step 7: Compute vd using (2.37). ' 
oi 
Step 8: Find {Hy} from (2.40). 


Step 9: Increase i by 1. 


Step 10: Calculate H using (2.42). 


Step 11: Goto Stepl. 


It may be aaded that the extra steps of the modified algo- 
rithm do not add significantly to the computation time of an iteration. 
Although faster converging error sequences have been observed when 
using the modified algorithm, an analytical proof of the fast conver- 
gence is not available at this time. However, this approach is useful 


in many optimization problems where the dimensionality is too high for 


the Fletcher-Powell algorithm [14] to be practical. 


i ita Seite 
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C. An Alternate Derivation 
It should be noted that the algorithm of Section IIA can be 
derived by applying a gradient [14] method of optimization to the error 


function 
“i-] -1 N-1 rm as rc 
een ye EME CH, I-A, 0 (2.50) 
mn kal ge0 k£ L 
Using (2.6), 
“i-1 
OH ky yPk yas 
ani! 
Pq 
and 
“i-1* 
one “ yk yas 
his? 
pq 
lead to ain] Afn4 
h = 
= n}) z M-1 N-1 Ay pk .-qg 
ai-] =| Dy = (1H, || TR) U V 
oh =0 g=0 H | 
Pq kg 
kg 
“i-1* 
kg BK G2 _ gi-l* pk ae 
eR gh ae Or By PE hl (2.51) 
His 


Equations (2.15), (2.16), (2.24), and (2.51) yield 


ai-1 
aK({h we 1) 
ai-] 
ah 
pq 


= -2MN(ho - f ye (2.52) 


pparlansendiisen anal 


tani 


att shana 
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If the classical gradient method were being used to decrease K({h*}), 


the next step would be to implement 


“i~] 
ni x pint, aK({h mn y ( 
pq pq fe] afiri ) p,q) eT, 
pq 
= 0, (p,q) €T, (2.53) 


where p is a negative number such that 


aK({ht 3) 
3p 


=0. (2.54) 


In the algorithm of Section IIA, (2.53) is implemented such that 


p=-1/(2MN) , (2.55) 


It can be shown that (2.53) and (2.55) do not give the result of (2.54). 
Therefore, the phase correction algorithm is related to, but not the 
same as, the classical gradient technique. 

Steiglitz, in [12], minimized an error function like that of 
(2.50) for one~dimensional recursive filters, using the Fletcher-Powell 
algorithm [14]. Unlike Steiglitz we did not use cascaded second order 
sections, but did use the FFT. We did not use the Fletcher-Powell 
algorithm because of storage considerations. In terms of real floating- 
point words, the storages of the original algorithm, the modified algo- 
rithm, and the Fletcher~Powell algorithm are 


Ss, = (M,+1)(N,+1) + 2MN (2.56) 


S, = 3(M,+ 1)(N +1) + 2MN 


Sy = ((My+1)0N +1))7/2 + 5 HM +1) 00 +1) + 2MN (2.58) 


respectively, where (M, + 1)(N, +1) is assumed to be even in (2.58). 


For M,=N,=7 and M=N= 32, 


= 4448, 


Here, the Fletcher-Powell algorithm requires almost double the storage 


of the modified algorithm. If M, N, M,. and N, are increased linearly, 


the amount of storage required for the i ete a algorithm in- 
creases much faster than the amount of storage required by the modified 
algorithm. For both algorithms FFT calculations take up most of the 
computation time. Therefore, for each iteration, the computation time 


is approximately proportional to MNlog, MN. 


D. Phase Sequences and Phase Error 
In this section a technique for choosing the initial phase 


sequence (o) is formulated and a method of comparing the phase 
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sequence {3 to a set of desired phase sequences is developed. 
The following notation and observations are needed in this section. 
The frequency response of the filter uh (2, .2,) is ex- 


pressed in amplitude-phase form as 
1 
Ha(exp(~jw ) exp(~jw 5) 
i 

= |HA (exp(-Ju 1) -exP(-Jw)) [exPGa (w) 2), (2.59) 
where lod (exp(-j,) ,exp(-jw,) )| is the amplitude response and 
aw, ,W,) is the phase response. Salient properties of the frequency 
response are: (i) The symmetry conditions for the amplitude and phase 


responses are 


[a (exp(-jw ,) ,exp(-jw 5) |= Ha(exp(jw )) expliw 4)) | (2.60) 


i i i 
a (w, -W,) =p (w, -w,)2n-0 (-w, +“) (2.61) 
where p(w, .W,) takes on integer values. Conditions (2.60) and (2.61) 
are necessary and sufficient for fn) to be a real sequence. 
(ii) Linear phase responses for digital filters may be expressed as 


(0, +09) = -aw,/2- pity /2+ §n/2+ vw, say) (2 .62) 


where @, 8, and § are non-negative integer constants, with § equal to 


0 or 1, and Vw, 15) is a function taking on values 0 or l. If §=1 


in (2.62), the filter is said to have negative symmetry and if § = 0, 
the filter has positive symmetry. 

(iii) The results of this section are derived for filters with positive 
symmetry but they are easily derived for negatively symmetric filters 
as well. 

If uh (2, .2,) has linear phase and positive symmetry, 


"(wy 1 Wy) may be expressed as 


1 i 
where K and L are non-n tive integer If i - oe 
; aoe iets m,N-1' “nd’ “Ok‘ 
and byl z are non-zero for some values of m, n, k, and t, the 


only possible values of K and L in (2.63) are M-1 and N-1l, respec- 
tively. 

(iv) It is specifically assumed in this section that T is described by 
(1.11), meaning that (2.3) may be rewritten as 


(2.64) 


if HA(Z 
te 


) is a linear phase filter with positive symmetry and if he 


Ge. ; 0’ 


i i 
mN,' Ao! and "Mt are non-zero for some values of m, n, k, 
and 4, the phase response Ow 5) of Hia(Z,.Z,) can be expressed 


as 


A" (ws +9) = -(M,/2)0, -(N ,/2)09 + ¥(w, +p) 0 (2.65) 


where Vw, 1W,) is a function taking on values 0 or l. 
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i i ai 
(v) The DFT that is used to get fH} from ie! and {H, ,} from 


(i) is of order (M, N). This means that W, and W, are restricted 


1 2 
by 
®) = (2n/M)k, 09<k <M-l, (2.66) 
w, = (2n/N)4, O54 eN-1, 2.67) 
and that 
Hy = Ha(2n/M)k, (27/N) 4), 
By 7 Ha ((2n/M)k, (2n/N) 4). 


Thus from (2.63), (2.66), and (2.67) a positively symmetric linear 


phase filter, Ha (2) .2,), has a phase sequence expressed as 


Og» = a ((2n/ Mk, (2/N) 2) 


= -(K/2)(2n/M)k - (L/2)(2n/N) £ 


+ ((2n/M)k, (2n/N) £) r. 
(2.68) 


Let 3 be a subset of Ty such that k, Lest if yiian/M)k, (2™/N)2) =0, 


gsi otherwise. Equation (2.68) can now be expressed as 


% 17 76/2) 2n/M)k ~(L/2) (2n/N) 4 (2.69) 


for &,2)es! and 
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%, = -(K/2)(2n/M)k-(L/2)(2n/N) 4+ (2.70) 
if (k,t)¢ st, It should be emphasized that a linear phase response 
implies a linear phase sequence; however a linear phase sequence 
such as (2.68) does not necessarily imply a linear phase response. 
For a filter Hi (2, .2,) to have the linear phase response of (2.63), 


the necessary and sufficient coefficient constraint is 


i 


i. « 
mn 


hy »>msKandnslL, 
-m,L-n 


«6 »m>Korn>L, (2.71) 


THEOREM 2: If the initial phase sequence (or) is described by 


(2.68), with (K,L) = (M,-N), every phase sequence (er) can also 
be described by (2.68) with (K,L) = (M,.N)). 


Proof: Let (K,L) = (M,.N,) and let U and V be as defined inChap- 
ter I. Note that 


exp(j %,) = ukM1/2 yN1/2 (2.72) 


if k, t)est and 


expla, = (2.73) 


if k,t)¢ st, For m=M 
(2.73) imply 


1 and ns Nie equations (2.4), (2.72), and 
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aj i a a) Z tH 

Bn 7 (1/MN) Mot NEB, lu k(m-M}/2),- ¢(n-N, /2) 
k=0 go “4 
(k, 4) est 


-(1/Mn) Met Nel IH, ju R(m-My /2),- e(n-N 1/2) 
(k, 2) 4S 
The fact that M,-m-M,/2 = -(m-M,/2) and N 


leads to 


1 


ai - M-1 N-1 k(m-M,/2),,¢(n-N ; /2) 
By mW nt O/MN) Og OE |B, |U Vv 
1 1 k=0 
kk, Des! 
~t Ne - -N,/2 
-(1/MN) M: 1 Nel is, oem M,/2),,4(n N,/2) 
= £0 
(k, 2) 4s! 
ai" .¢! 
mn mn 


for m $s M, and ns N)- Using (2.68), (2.71), and (2.75) 


af 
é.,° -(M, /2) (2n/M)k (N /2)(2n/N) & 
+ 9) ((2n/M)k, (2n/N) On 
Equations (2.21) and (2.76) imply 


Og p= (My /2) (2n/MIK= ON, /2) (2n/N) 4 


+ 7 4((2n/M)k, (2n/N) £)n 


(2.74) 


“n-N,/2 = ~(n-N,/2) 


(2.75) 


(2.76) 


(2.77) 
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Replacing 94( (2n/M)k, (27/N)2) by yt! an/M)k, (20/N)2) in @.77}, 
tor} can be expressed by (2.68) if (or) can be expressed by 
(2.68), when (K,L) = (M,,N,). Since {80} is described by (2.68) 
with (K,L) = (M,.N,). the theorem holds by induction. 


Q.E.D. 


Corollary: If (e] is described by (2.68) and (K,L) = (M,.N,), the 


error sequence W converges ina finite number of iterations. 


Proof: Assuming (K,L) = (M,/N,) and (2.68) holds for i=0, theorem 
2 implies (2.68) holds for every i. The function vi((an/M)k, (27/N)2) 
takes on one of two values for every one of the MN ordered pairs (k,?) 
in Tp ‘ Therefore there are only gus jens candidates for the func- 
tion VY ((2T/M)k, (2™/N)t) and fewer than 2 when (2.19) holds. 
Since W is strictly decreasing before convergence and since different 
error values - correspond to different phase sequences (or, there 


are fewer than il iterations. 


Q.E.D. 


The corollary to theorem 2 implies that if (e? ,3 is properly 
chosen the phase correction algorithm may be used to design linear 
phase response filters. However, in designing such filters, numerical 
problems in the FFT implementation of the DFT cause the symmetry 
condition of (2.71) to be violated. In other words 


i 
Rian ~ Dy an <4 “gyt SE ees 


i 
where ane is an error sequence with non-zero members. The result 


+ 
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i 
is that the phase sequences {e ) may not be described by (2.68) and 
W does not converge ina finite number of iterations. This can be re- 
medied by implementing the arithmetic replacement statement 


= 5(ht eh ie 


In general the limit og of the sequence W depends upon 
the initial phase sequence ies as illustrated in the section on 
numerical results. In other words, the algorithm converges to a local 
minimum rather than to a global minimum. It may be possible to gua- 
rantee convergence to a global minimum in the 1-D case by making 
use Of minimum phase versions of the filter coefficient polynomials. 
However, polynomials of two or more dimensions generally do not have 
minimum phase versions and such a procedure cannot be extended to 
two or more dimensions. At present, the only way to insure that this 
error I is small is to pick (o) such that :" is small. If ae is 


small, the fact that W is non-increasing precludes the possibility 


aoc 
that I is large. 


Let 
0 = - 
Beg = 7 (20/M) K/2)k-(2n/N)(L/2) 44Y, (2.78) 
0 for K and L even, 
u(k-M/2) for K odd and L even, 
Vue = 


u(t-N/2) for K even and Lodd, 


u(k-M/2) + u(-N/2) for K and L odd, (2.79) 


0 forx <0, 
u(x) = 


l forx 0, 


where K and L are non-negative integers which are to be chosen such 


| that yo is minimized. Typically, ‘hi is minimized for some (K,L) 


close to (M,,N,). However, if (K,L) = (M,,N)), there are a finite 


number of iterations, using the corollary to theorem 2, and the limit 


of W is not much smaller than ‘ee For this reason one should not 


use (K,L) = (M,N ) in (2.78) unless a linear phase response filter 


1 


is mandatory. 


Theorem 3: If (e?,) is specified by (2.78), then 


Wo m- . 
neo L2M-KNoyqel2n-Llony 


Fat lH,, , [exe i7Yy 4) , (2.81) 


where lH, =0 fork 2M or t2N, and where the DFT is for 


sequences of length (2M,2N). 


Proof: From (2.9) and (2.17), 


m=0 n=0 


L 
Equation (2.78) leads to exp(i®),) = oy i/2 exp (iT, ») 


and therefore, using U = exp(-j27/2M) and V = exp(-j27/2N), 


-k (m-K/2), g(n-L/2) 


exp(jTY, ») 
2M- = cv-k(2m-K) ~~ g(2n- 
= (a/amany) 7M Nig | oem k yo s2o-L) amy, |) 
k-0 fo“ - 
(2.83) 


ea T[2m-K), .,/[2n-L., 


toh (i ; 
On substituting *9(2m-Kl, for hon in (2.82) one obtains 


ui 2n-Lo ny 


the result. 


OLE. D.. 


After using the IDFT to find fan! One can easily find from (2.80) the 
values of K and L for which 199 is minimized. 

With the exception of the special case of theorem 2 and its 
corollary the phase sequences (9) produced by the phase correction 
algorithm are in general nonlinear. Linear phase sequences are usually 
desirable and it is therefore profitable to compare the actual phase se-~ 
quences (e3 to linear phase sequences. Let the linear phase se- 
quences, which we shall compare the actual phase sequences to, be 
described by 


d 
9 


aa" 6m—(2n/M) (K/2)k-(2n/N)(L/2) £47, » (2.84) 


and (2.79) where K, L, and 6 are non-negative integer variables and 
6 is 0 or 1. A convenient measure for the error between {o,,3 and 


the linear phase sequences of) is 


M-1 N-1 2 d South 
r EW, | Jexpge, )-expe, 1). 2.85) 


J,=min  ( 
6,K,L k=0 


i 


ars 
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Two advantages in using J as a measure of phase error are: 

(i) phase differences of 2 are equivalent to zero phase differences 
and 

(ii) phase errors in pass regions are weighted more heavily than phase 
errors in transition and reject regions. 


Suppose that fer some positive A smaller than 7, 


d i 
lm 47 % 4 =A, (k, g) eT) (2 . 86) 


Equation (2.86) implies that the error between (er 1} and a linear 
phase sequence equals 4 at each of the MN DFT points in the fre- 


quency domain, and that 


i d 
Op Pay? O° (k, 2) eT): (2.87) 


If (2.86) holds and (2.87) is substituted into (2.85), the fact that 


2 
JexpU ay )-expU (of + »)) | = 2(1-cos(4)) implies that 
- 2 
l= "3 zr |H, | 2(1-cos(4)) 


M-1 N-1 2 
2(1-cos(s)) x 5 |Hy | 


k=0 g=0 


2MNE(1-cos(,)). (2.88) 


For a given value of A, if 


J, < 2MNE(1-cos(a) (2.89) 
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then 
d i 
I 47 yl #4 (2.90) 


will probably not be true for all (k,%) in To: but it will tend to be true 
in the pass regions, since the error in the pass regions contributes 
heavily to the value of Jie 

Given that (2.90) is desired for some positive 4 <7, J, 
can be calculated every few iterations and the algorithm can be halted 
when (2.89) is not satisfied. If the initial phase sequence is a linear 
one such as given in (2.78), (2.89) will hold at least for i=0. 

It is assumed that the specific values of 6, K, and L in 
(2.84) are not of critical importance. If the values are needed, they 
are easily determined by using the following theorem. 


Theorem 4: For given values of E and Lame 


J, = 2MN(E-max(/R_1)) 


i 
lH, for K and L even, 


[nt UK? expuny, ,) for K odd and L even, 


4/2 


i 
ben ee By exp(jTY, ») for K evenand L odd, 


/ 
k/2 ft, exp(imy, ,) for K and L odd. 


i 
Ran <> Hy, | HU 


m 


j Proof: We prove the theorem for K and L even. The proofs for the 


i other cases are similar. From (2.84), (2.85), and (2.9) 


| a 1 2 M-1 N-1 
=| |. Le ee JH, |-2max § Ns IH, Hy, exP- i9p,)) 
k=0 es 6,K,L  k=0 ae 


ai a 
= 2(MNE-max ((-1)§ M5) NS? iw, HK, 2 sin iced 
6. K,L 6 a 
Replacing (K/2, L/2) by (m,n) leads to 
-km vy) 


J, = 2MN(E-max (|(1/( MN)) = SH, AL 
m,n k=0 ¢0 


Q.E.D. 


In theorem 4 we show how to measure phase error ina 
filter. In the following chapter we show how to measure the effects 


of this phase error on filtered pictures. 


III. THE MEASUREMENT OF PHASE DISTORTION 


A. Derivation of the Measure 
In this section we develop a measure of phase distortion, 
and give a few of its properties. Let {pnt be an output picture, with 


M rows and N columns, satisfying 


p <> Pie (3.1) 


} 


and having a Fourier transform P(wW,, W,) as described in (1.17). {Py 
is a sampling of one period of P(w) .W,) at MN points. {P,! satis- 


fies the radial symmetry condition 


P (3.2) 


ee PIM-K) [N-4, 


The output picture eee which has the Fourier transform Pw, w,) 
described in (1.19), satisfies 


Aes a ee (3.3) 

Bgl = egh (3.4) 

Pe P, pexP(im, ») ; (3.5) 
and 

m4 = m((2™/M)k, (2T/N)L) , (3.6) 
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where m(w, ,W,) is given in (1.20). 
: A measure of phase distortion in {pn due to the effects 


of m(W),W,), is 


M-1 N-1l 2 
Ak ae re aon neo ono me 
D(p,p) my M-l N-1 

2 


» x 
‘ m=0 n=0 Pian 


(3.7) 


which is the normalized squared Euclidean distance between pat 
and {Pont The measure D(p,p) can be expressed in the frequency 


domain as 


M-1 N-1l 


i) 4 
es: Z 
Sc Re ee ee ; 
D(p,p) > M-1 N-1 3; ) 
2 
5 
2 2) By |Pye 


Filtering shifts the output with respect to the input. The fact that a 
shifted picture is recognizably the same as one which has not been 
shifted means that we must remove the effects of shifts due to filter- 
ing from D(p,p). The phase response of a filter which linearly shifts 
the output with respect to the input is 


2(w,,@,) = ~ (K/2)w, - (L/2)w, (3.9) 


where K and L are non-negative integers. The phase sequence {¢, ats 
derived from (3.9), is 


G4 = ~ (K/2)(27/M)k - (L/2)(27/N)L + Why (3.10) 


where Y,, is given in (2.79). 


kt 


From (3.10) and the definitions of U and V, 


kK/2 yl/ 2 


exp (iP, ») =U exp(jTY,,) . (3.11) 


kt 
In addition to its shifting effect, filtering can cause the output picture 
to be close to or equal to the negative of the input picture. We can 
remove this effect from D(p,p) by using a variable § which equals +1. 
Based upon these observations our measure of phase distortion is 


M-1 N-1 x 
i z x © - 
min | Py gexPliP, ») oPe 


a 6, K,L k=0 2=0 
D(p,p) eee eT oe ee 


(? 


This measure D(p,p) has the property 
0 = D(p,p) 1 


Two advantages to using D(p,p) as a measure of phase distortion are: 
(i) phase differences of 27 are equivalent to zero phase 
differences and 


(ii) the amplitude sequence weights the phase differences. 


B. Evaluation and Application of D(p,p) 
In this section we show how to evaluate D(p,p) and how to 


use it in judging filter phase responses. 
e 


{PB ,}, M, and N, 


Theorem 5: Given {P Khe! 


WEE 


D(p,p) = 1- (maxlb  )/E 


RO RRR Ponce onerns mm nn 
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where 
M-1 N-1 2. 
=(1/MN) = = |p|” , (3.13) 
k=0 4=0 
boo Pe Pe for K and L even, 
k/2, 4/2 
By ss << PsP 0 V exp(jTY, ») for K and L odd, 
bo <> P* Pp U2 cup (inv ) for K odd and L even, and 
kt kt kt 
bn <> p* vi’? (j7tY, ,) for K even and L odd. (3.14) 
eee Pus¥ expQry,, . 


Proof: The proof is for K and L even. The proofs for K and/or L odd 
follow similarly. Equation (3.12) yields 
M-1 N-1 


2 a 
min = z IP, | + eel ~ 26PR) P, pexp(- J, ») 
Pe Fe 10 eee aris SUL eens a ei NS 
D(p, Pp) = oe 1 2 1 5 
@ 8 ak, they! 


which can be simplified, using (3.11) and (3.13), to give 


M-1 N-1l 
E- max §((1/MN) 2 py ee = ye 
6, K,L k=0__4=0 


E 


D(p,p) = 


Replacing K/2 by m and L/2 by n and using (3.14) for K and L even 


' we have 


D(p,b) = 1 - (max, (6b) )/E 
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The variable § can be eliminated by an absolute value, giving 


Dip.) = 1 - (maxlb,|)/E 


Q.E.D. 


Given an example picture {pa and several sequences, 
{m, ,3, of phase distortion numbers, equations (3.3) - (3.5) can be 
used to generate degraded pictures ae and the corres ponding values 
of D can be calculated using theorem 5. With the values of D and the 
degraded pictures, we can determine how large D can get before the 
picture {Bn is unrecognizable. We denote the largest allowable 
value of D(p,p) as Dia 

We have shown how to calculate D(p,p) when the output 
picture frequency characteristics {PL Pe and - are given. In the 
following we show how to calculate D(p,p) when the frequency charac- 
teristics of the input picture and filter are given. The DFT versions of 
(1.17) and (1.19) are 


= ad 
Pee = Fyy| Hyp lexptiey,) (3.15) 
and 


Pp Fy! Hy, lexpdie, ») (3.16) 


where 


27 H((27/M)k, 27/N)t) , 


=" 


| 


00 = OTC (2M/M)k, @M/N)L) , 


and 
) “= 8((2™/M)k, (27/N)2) 


k 


Notice that the factor P*,P occurs in each of the expressions of 


KL kt 
(3.14). Using (3.13), (3.15) and (3.16), 
i M-1 N-1 4 ; 
E=(1/MN) = £ le) Can (3.17) 
k=0 t=0 
and 
ah aia 2 2 d 
Bey = lr, I ln, ,| exp(j(®, »-8 »)) (3.18) 


It is interesting to note that the amplitude of {Ft but not its phase, 

is used in (3.17) and (3.18). D(p,p) may be calculated using (3.17), 
d 

(3.18), and theorem 5, given only Ce, | is flu, ,| }, WEE and 

{oe 


sponse, we may calculate D(p,p) for the filter and several input pic- 


et - To find out whether or not a filter has an acceptable phase re- 


tures or test patterns. If D(p,p) tends to be greater than Da the filter 
is not useful. This testing procedure uses much computation time and 
storage if M and N are large, 

Equations (3.5), (3.15) and (3.16) correspond to circular 
convolution and theorem 5 requires the evaluation of a circular correla- 
tion, Also, the frequency responses exp(im(w,,W,)), Hw, ,W,) and 
H(w) ,,) generally correspond to an infinite number of filter coeffi- 


cients, implying that M and N should be infinite. However, our main 
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purpose is to find and use D(p,p) rather than to filter pictures, so 


circular convolution and circular correlation are acceptable. 


C. Estimation of D(p,p) 
When sampling H(w,,W,) to produce {HL t, M and N 
may equal 128, for example, and be large enough for the representa- 


tion to be accurate. However, wren [F(w w,)| and PW), W,) are 


1’ 


sampled to produce SAE: and {P. M and N are typically equal 


a3 
to the picture's dimensions, which may be much larger than 128 by 128. 
If the input spectrum [Pw ,w,)| is modeled by a smooth function 
[F(w,.,)| , which is accurately represented by a small number of 
samples, we can greatly reduce M and N and estimate D(p,p), 
thereby saving storage and computation time. 

The Fourier transform of a one-dimensional pulse is a sinc 
function which dies off as 1/W and has a value of 1 for W=0. Simi- 
larly [15] the two-dimensional Fourier transform of a disk shaped pulse 


dies off approximately as 1/R where 


2,1/2 


9) (3.19) 


R= (w* +W 
1 
A simple model of the input picture amplitude response [Fw .w,)| is 


|F(w, ,w,)| = 1/(1 + aR) (3.20) 


which dies out approximately as 1/R. The constant a may be varied 
until the approximation is accurate. The digital version of IFiw, ,,)| 
is Clea I, which is defined by 


Ir. ,L=1F@,.%,)1 , 


semen 
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(21/M)k, k*M/2 


(2m/M)k-27, k>M/2 , 


(2m/N)t, tS N/2 
(2M/N)t-27, L%>N/2 . 


Replacing cE pt } by iF, in (3.17) and (3.18), we may use theorem 
5 to calculate D(p,p), an estimate of D(p,p). In general [Fw,.,)| is 
a much smoother function than IF, .w,)| , allowing us to reduce M 
and N, save storage, and save computation time. 

The estimation technique of this section greatly simplifies 
the testing procedure for filter phase responses of section IIIB. Assume 
that [Fw, ,w,)| is a model for all input picture amplitude spectrums 
|F(w, ,,)| of interest to us. To test the phase response of a filter we 
can compute D(p,p) one time and RomEetS it to ae When a filter is 
designed by the method of Chapter II D(p,p) may be calculated at each 


iteration. When D(p, p) becomes larger than Di! the design algorithm 


may be terminated. 


IV. NUMERICAL RESULTS 


A. Filter Design Examples 

The numerical implementation of algorithms described in 
Chapter II, together with the design of filters and a comparison of their 
frequency responses is presented in this section. The computation was 
performed on the CDC-6600 at The University of Texas at Austin. The 
pictures are printed in 16 grey levels with darker levels representing 
more positive quantities. For every grey level picture the range be- 
tween the maximum and minimum intensities, of the function being dis- 
played, was linearly divided into the 16 levels. Except for the linear 
phase response filters, every filter was designed in 100 iterations, 


with i = 99 and I_=0. 
max E 
Example 1 The desired amplitude response is a ring described by 


CNY las 
2) : 


2 
= < 


(4.1) 


-02 otherwise, 


1' Ny are 32, 


32, 7, and 7, respectively. The amplitude sequence (lH, ,13 is shown 


where jw, | s 7 and |w, | stm, The values of M, N, M 


in Figure 3. Using the initial phase sequence of (2.78) with K and L 
determined from theorem 3 to be 8, Filter 1.1 was designed by the algo- 
rithm of section IIA and Filter 1.2 was designed by the modified algorithm 
of section IIB. The amplitude and phase responses of Filter 1.2 are shown 
in Figures 4 and 5, respectively. The responses of Filter 1.1 are almost 


identical to those of Filter 1.2, and are not given. The coefficients of 


Filters 1.1 and 1.2 are presented in Tables 1 and 2. The initial and 
final error numbers for Filters 1.1 and 1.2 are normalized using the 


number E and given by 


00 rey = (°° /g) = 213185 , 


(199 99 7p) ~ 10368 , 


(799-99 75) ~ 102859 . (4.4) 


One gets an idea of the rate of convergence of ae and lag as i799 
from the plot of log (tye) - ,1028) and log (tye) - .1028) versus 
i, given in Figure 6. As the figure shows, the modified algorithm's 
error sequence Ww is nearer convergence after 40 iterations than the 
Original algorithm's error sequence W after 100 iterations. The coef- 
ficients of Filters 1.1 and 1.2, in Tables 1 and 2, are significantly 
different. Filter 1.3, also with the desired amplitude sequence of 
Figure 3, had the linear initial phase sequence of (2.78) with (K,L) = 
(M,.N,) = (7,7). The filter was designed using the original algorithm 
and its amplitude and phase responses are given in Figures 7 and 8. 
The coefficients are given in Table 3. The error sequence W, for this 
filter, converged at the fourth iteration, which is in agreement with the 


corollary to theorem 2. The normalized error numbers for Filter 1.3 are 


moe = .1374 (4.5) 
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(33g) ~ 1302. (4.6) 


As predicted by theorem 2, the phase response of Figure 8 is linear. 
The error numbers of (4.5) and (4.6) show that H3 (2, ,2,) is not sig- 
nificantly better than Ae (2) .2,). The final filter is usually a small 
improvement over the original filter when (2.78) holds and (K, L) = 
(M, N))- 

Comparing Figures 4 and 7, one can see that Filter 1.2, 
with the nonlinear phase, has a significantly smaller magnitude of 
ripple in the reject regions. Whereas Filter 1.3 has four distinct 
ripples in its passband, the passband ripple of Filter 1.2 is not de- 
tectable in the figure. The passband rejectband boundaries are dis- 
tinctly diamond shaped for Filter 1.3 but are more nearly circular for 
Filter 1.2. Comparing Figures 5 and 8, one can see that the phase 
response of Filter 1.2 is not wild in the passband, but fairly linear. 

Filter 1.4 was obtained by windowing Filter 1.3. Each 
coefficient oa of Filter 1.3 was replaced by he Won in Filter 1.4, 


where 
to = i= dee ; 
Fenn = Em = M,/2)° + (n= N,/2)7) 7? 
R= [(M,/2)” + (N,/2)°J 17? _, 

and 


Doe 
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The window in this example is a generalization of a window found in 
[8]. The amplitude response of Filter 1.4 is shown in Figure 9. Com- 
paring Figures 7 and 9, it is apparent that although windowing reduces 
amplitude response ripples it also decreases the steepness of the re- 


sponse in transition regions. 


Example 2 For this example the desired amplitude response is 
= s 9 
A(w,,@,) =1, lw, | + lw, $2.1 
=0, otherwise (4.7) 


where Jw, | <7 and Jw, | st, This corresponds to a low-pass filter 
with a diamond shaped passband. Both filters of this example were 
designed by the modified algorithm with M, N, M, and Ny equal to 
64, 64, 7, and 7, respectively. The desired amplitude sequence for 
the filters is shown in Figure 10. 

The initial phase sequence of Filter 2.1 was (2.78) with 
(K, L) = (8,6), using theorem 3. The amplitude and phase responses 
of Filter 2.1 are given in Figures 11 and 12 and the coefficients are 


presented in Table 4. The normalized error numbers for Filter 2.1 are 


(79° /e) © .08691 (4.8) 
and 


799.99 


( /E) * .08198 . (4.9) 


For Filter 2.2 the initial phase sequence was that of (2.78) 


with (K,L) = (63,63). The amplitude and phase responses are displayed 
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in Figures 13 and 14 and the coefficients are presented in Table 5. 
Plots of (la | } and (1a? | } for this filter are given in Figures 15 


and 16. The normalized error numbers for the filter are 


(°° /e) * .99852 (4.10) 


(799-99 7p) = 10098 (4.11) 

It is apparent from (4.10) and (4.11) that the final filter is 
greatly improved over the initial filter. Figures 15 and 16 show that 
most of the coefficient energy is concentrated in those points corres- 
ponding to the set T. However, (4.9), (4.11), Figure 11 and Figure 13 
indicate that Filter 2.2 has a much larger error number and a worse am- 
plitude response than Filter 2.1. Also, ie for Filter 2.1 is smaller 
than I~’ for Filter 2.2. This illustrates the value of choosing a 
good initial phase sequence. 

The values of M=N=32 in example 1 and M=N=64 in 
example 2 represent fairly coarse frequency sampling. For the exam- 
ples we have chosen, however, larger values of M and N make no 


substantial] difference. 


B. Phase Distortion Examples 
In this section we distort the phase of a test pattern by the 
method of (3.3 - 3.5) and by filtering and we then measure D for the 


degraded images. The actual values of D are compared to estimates 


obtained through the modeling procedure of section IIIC. 


{ERATE 
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Example 3 Using (3.3 - 3.5) we distorted the phase of the spoke test 
pattern, {pon of Figure 17 with three zero mean Gaussian noise se- 
quences {m, zs . The sequences differ only in their standard deviations 
(i.e., the sequences are multiples of each other). The resulting pic- 
tures {Pa for standard deviations of .13, 375 and 1.0 are shown in 
Figures 18 through 20 respectively. The corresponding values of D(p,p) 
are: .0096, .0787 and .32. 

It is apparent from the figures that as a picture's phase 1s 
distorted its energy is redistributed in the spatial domain. This is a 


type of smearing. 


Example 4 As a second example we distorted the phase of {pint in 


Figure 17 using the deterministic function 
m(W.,W_) = 9(w* + we) 1/2 (w +w)., (4.12) 
i 2 : 1 2 1 2 


The resulting picture {pnt , shown in Figure 21, corresponds to 
D(p,p) = .0779 . (4.13) 


It is apparent that although the deterministic distortion of Figure 21 
differs from the random distortion of Figure 19, the figures are equally 
recognizable as spoke patterns and the values of D for the two figures 


are about the same. 


Example 5 Asa third example we filtered the test pattern of Figure 17 
using Filters 1.2 and 1.3. The desired amplitude response for both 
filters is given in (4.1) and the desired phase response is that of (3.9) 
with K=L=7, 
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For Filter 1.2 the actual value of D, with M=N = 128, is 
D(p,p) = .078 (4.14) 


and the output picture is shown in Figure 22. For Filter 1.3 the actual 
value of D, with M=N = 128 is 


D(p,b) = .00974 (4.15) 

and the output picture is shown in Figure 23. 
From (4.14), Filter 1.2 produces about the same amount of phase dis- 
tortion as that seen in Figure 19. From Figure 19 one can correctly 
infer that the phase distcrtion seen in Figure 22 is moderate. If op- 
posing spokes in Figure 22 are compared it is noticed that, unlike the 
spoke pattern of Figure 17, this filtered spoke pattern lacks radial 
symmetry. From (4.15), Filter 1.3 produces almost the amount of phase 
distortion seen in Figure 18. Although Figure 18 reflects a noticeable 
amount of phase distortion it is not easily detectable because Filter 1.3 
has linear phase and Figure 23 has radial symmetry. 

For both filters, the values of D(p,p) were estimated using 
the method of section IIIC with M = N = 32 and a = 32. In Figures 24 


and 25 respectively, Clrglt and our model {lp |3 are shown for 
M=N= 128. a is picked to be 32 for the model so that the dark re- 


gions in Figures 24 and 25 are approximately the same in area. For 


Filter 1.2 our estimate of D is 
D(p,p) = .07334 


while for Filter 1.3 our estimate of D is 


D(p,p) = .00926 . (4.17) 


Comparing (4.14) to (4.16) and (4.15) to (4.17) we see that the estima- 


tion procedure of section IIIC works well for our examples. 


V. CONCLUSION 


By allowing the desired phase sequence of a filter to vary, 
while the desired amplitude sequence is held constant, one can signi- 
ficantly improve the filter's amplitude sequence. The phase correction 
algorithm always converges and its convergence can be considerably 
speeded. Because of the simplicity of this design method, it is rela- 
tively fast, with good filters being designed in less than 30 sec. for 
M = N = 32. Through the use of theorem 3 good initial phase sequences 
may be picked. An easily applied DFT correlation theorem makes it 
simple to measure the phase error between the actual filter and a linear 
phase version of it. 

As our examples show, it is possible to measure the phase 
distortion, due to filtering, in digital pictures. Using the measure, 
the acceptability of filter phase responses can be determined. Using 
a simple model of input picture amplitude responses, the measure D 


may be estimated. 
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Figure 1. Block Diagram for the Original Algorithm 
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Figure 3, {]H, 1) For Filters 1,1, 1.2, and 1.3, 
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Figure 4, Amplitude Response of Filter 1.2. 
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Figure 5. Phase Response of Filter 1.2. 
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Figure 7. Amplitude Response of Filter 1.3. 
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Fig. 9. Amplitude Response of Filter 1.4. 
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Figure 10. {(H, ,[} For Filters 2.1 and 2.2. 
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Figure 11. Amplitude Response of Filter 2.1. 
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Figure 12. Phese Response of Filter 2.1. 
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Figure ]°. Amplitude Response of Filter 2.2. 
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Figure 14. Phase Response of Filter 2.2. 
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Figure 21. (Bon! for D(p,p) = .0779 
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Figure 22. Output Picture for Filter 1.2 
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